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By Philip R. Nachtsheim 
Ames Research Center 


SUMMARY 


Oblique wave formation in the melt layer of a body entering the 
atmosphere is examined. The three-dimensional disturbance equations of the 
liquid including the effects of viscosity are formulated. The disturbance 
motion of the gas is taken into account neglecting viscosity. The Tollmien- 
Schlichting type of instability in the liquid is excluded and only those modes 
are considered for which the phase velocity of the disturbance is greater than 
the velocity of the basic liquid flow at the gas- liquid interface. Numerical 
and approximate analytic results are presented for waves in the liquid with 
supersonic flow over their crests. A comparison with the numerical solutions 
shows that the approximate analytic solutions are accurate for solving the 
eigenvalue problem when the Reynolds number of the liquid film is small, that 
is, for a highly viscous melt layer. The approximate solutions are obtained 
in closed form. It is demonstrated that the mechanism of liquid wave genera- 
tion is supersonic wave drag. The relation between the wavelength of the 
disturbance and the other parameters involved in the eigenvalue problem is 
determined and is presented for a wide range of these parameters. 


INTRODUCTION 


The theory of hydrodynamic stability has been employed by several authors 
in analyzing the formation of waves in the melt layer of a body entering the 
atmosphere. 

Feldman (ref. 1) considered a liquid shear flow with wave crests parallel 
to the stream direction (streamwise grooves) but neglected the disturbance 
motion of the gas. Miles (ref. 2) considered a liquid shear flow with wave 
crests normal to the stream direction, and although he did not include the 
disturbance motion of the gas in calculating the stability of the liquid film, 
he did indicate that the disturbance motion of the gas cannot be relegated to 
a subsidiary role when the phase velocity of the disturbance wave is greater 
than the velocity of the liquid at the gas-liquid interface. Chang and 
Russell (ref. 3) considered wave crests normal to the stream direction and 
they included the disturbance motion of the gas in their analysis. However, 
the undisturbed liquid configuration considered by them was the classic Kelvin- 
Helmholtz type, as treated by Lamb (ref. 4), namely, an infinitely deep, 
initially quiescent liquid suddenly subjected to a disturbance, periodic in 
time . 



The waves considered in this paper are the type observed, for instance, 
in reference 5, The crests of the waves are oblique to the stream direction 
(fig. 1); the problem considered is the three-dimensional disturbance of a 

shear flow of a thin liquid film adjacent 
to a high-speed gas stream. The distur- 
bance motion of the gas is accounted for in 
a manner similar to that employed by Chang 
and Russell (ref. 3). 

The problem of oblique (three- 
dimensional) wave formation cannot, in gen- 
eral, be transformed to a two-dimensional 
problem by employing Squires* theorem 
(ref. 6). For example, in the case of a 
supersonic gas stream, the gas flow over 
the crests of the waves in the liquid film 
may be either subsonic or supersonic, 
depending on the direction of wave propaga- 
tion. In subsonic flow, the disturbance 
pressure of the gas is in phase with the 
interface whereas in supersonic flow it is 
in phase with the slope of the interface. 
Clearly, the stability problem is different 
for the two cases and therefore requires a 
three-dimensional treatment. The primary 
interest in the present investigation is 
the coupling between the disturbance motion 
of the gas and the stability of the liquid. 
Specifically excluded are modes of energy 
transfer as treated by Miles (ref. 2), that 
transfer energy from the mean motion of the liquid to its disturbance motion. 
As Miles pointed out this type of energy transfer occurs in the Tollmien- 
Schlichting type of instability where the phase velocity of the disturbance 
wave is less than the velocity of the liquid at the interface. Hence, to 
study the effect of the disturbances in the gas on the stability of the liquid 
film, consideration will be given to those possible modes mentioned by Miles 
(ref. 2) for which the phase velocity of the wave is greater than or equal to 
the velocity of the liquid at the interface. Of course, considering these 
"fast" waves in the liquid raises the question of the Tol lmien-Schlichting 
type of instability in the gas, but this question can be ignored since in the 
majority of practical applications the gas flow is likely to be turbulent. 

In the following sections the three-dimensional disturbance equations are 
formulated for the liquid and the gas, and for the boundary conditions. The 
equations are then made dimensionless . 

A numerical method of solving the eigenvalue problem and an approximate 
analytic solution are presented. The approximate solutions are compared with 
the numerical solutions. Then neutral curves resulting from the numerical 
solution are presented for a selected range of the pertinent parameters. 
Finally, the approximate solutions are presented in dimensional form for a 
wide range of parameters for the small Reynolds numbers of the liquid film. 


External flow 

♦ 



Figure 1.- Two sets of oblique waves. 
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FORMULATION OF DISTURBANCE EQUATIONS 


Consider the stratified motion of a thin film of liquid over a solid 
surface adjacent to a high-speed gas stream. For analyzing the stability of 
this flow configuration, a perturbation is introduced into the describing equa- 
tions for the liquid and the gas. The basic liquid flow is considered to be a 
shear flow, and viscous effects are included in formulating the disturbance 
equations of the liquid. The effect of viscosity is neglected in the distur- 
bance motion of the gas. The only effect of the gas considered in the stabil- 
ity analysis is the pressure variation of the gas at the disturbed interface. 

Gradients in the basic flow of the gas are neglected except at the 
interface where it is recognized that the shear in the gas balances the shear 
in the basic flow of the liquid. Away from the interface, the basic proper- 
ties of the gas flow that influence the pressure disturbance are characterized 
by parameters regarded as constants throughout the gas boundary layer. The 
purpose of this investigation is to determine how those parameters affect the 
stability of the thin liquid film. 

Cartesian coordinates are introduced so that the interface coincides with 
the (xj, x 2 ) plane. The film of liquid flows over a solid wall, located at 
X3 = -h, so that the X3 axis points away from the wall. Cartesian tensor 
notation is used. The subscripts a and 3 range from 1 and 2 (two mutually 
perpendicular directions in the plane of the interface), and repeated sub- 
scripts imply a tensor summation over this range. The basic flow is taken to 
be in the xj direction; however, a three-dimensional disturbance is consid- 
ered. Furthermore, the wave amplitude is assumed to be sufficiently small that 
the principle of superposition applies. Therefore, the two families of wave 
crests shown in figure 1 may be treated independently. 


DISTURBANCE EQUATIONS - LIQUID 


Consider a three-dimensional disturbance of the basic flow whereby only 
first-order terms are retained: 


U a ( x l, x 2> x 3> O = U a (x 3 ) 

+ u a (x 3 ) 

E 

1 


u 3 ( X 1 > x 2> x 3> t) = 

u 3 (x 3 ) 

E 

> 

(i) 

P ( x l, x 2> x 3 > t) = P(X 3 ) 

+ p( x 3) 

E 




where 


E = exp i(kgxg - mt) 

The wave numbers k 1? k 2 are taken to be real. The wave normal has the 
direction numbers 

(ki, k 2 , 0) 
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The wave length of the disturbance in the direction of the wave normal is 

given by 2n/k where k = The constant to is taken to be complex. The 
real part of to is equal to the angular frequency of the disturbance, and the 
imaginary part is equal to the time amplification factor of the disturbance. 


The appropriate form of the disturbance equations is derived by 
substituting the form of the disturbance introduced above into the equations 
of motion of a viscous liquid 


3u 3u 

a — a 

3t + U $ 3xg 


^3 


3u 

a 

3x 3 


IiL + 

P 9x a 



3u 3 __ 3u 3 3u 3 

W~ + U 3 + n 3 93 ^ 


I 9£_ 

p 8 x 3 




> ( 2 ) 


j 


5u 3 

9x 3 


0 


and retaining first-order disturbance quantities. (Symbols are defined in 
appendix A.) The basic flow satisfies 


u l( x 3) = U x (l + x 3 /h) 


U 2 (x 3 ) = 0 


-h x 3 <_ 0 


P C x 3 ) = Pg - Pg x 3 


(3) 


where Pg is the pressure at the undisturbed interface. It is considered 
sufficiently general to allow the gravity vector to be normal to the (x^, x 2 ) 
plane. The disturbance equations are as follows: 


dU a 

i(kgUg - 0 ))u a + U 3 



+ V 




ifkgUg - w)u 3 


I + 

P dx 3 


v 



(4) 


iCkgUg) 



0 


J 
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In obtaining the disturbance equations, the fact that the basic flow 
satisfies equations (3) has been utilized. When u a (a = 1, 2) is eliminated 
from the foregoing equations the disturbance pressure p is expressed in 
terms of u 3 



where 



a 


is the velocity component of the basic flow in the direction of the wave 
normal. Eliminating p from the resulting equations yields the usual Orr- 
Sommerfeld equation for the liquid film 


<!!hi . 2k 2 „ 

dx 3 dx i 


k 4 u 3 





( 6 ) 


For a linear profile, d 2 U/dx 2 


0 . 


DISTURBANCE EQUATIONS - GAS 


gas. 


The effect of viscosity is neglected in the disturbance motion of the 
The following three-dimensional disturbance is considered: 

U a (Xj, X 2 , x 3 , t) = U a + u a (x 3 ) E 
usCxj, x 2 , x 3 , t) = u 3 (x 3 ) E 

> 


(7) 


P = P g + P(x 3 ) E 

p- = p g + p(x 3 ) E , 


where 


E = exp i(kgxg - wt) 
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and 


P = 




Y 


P 


oo 


Since viscosity is neglected in the disturbance motion of the gas, the 
appropriate form of the disturbance equations may be obtained by substituting 
the form of the disturbance introduced above into the equations of motion of 
a compressible, inviscid gas. 


^ ^ u a . — 8u a 

at + u s axg + u 3 ax 3 


a 2 3p 


P 3x a 


A 


9u 3 _ 9u 3 _ 3u 3 _ a 2 3p~ 

at + 3xg + U 3 3x 3 ax 3 


Bp _ 3p 
3l +U ^ = 



J 


( 8 ) 


and retaining first-order disturbance quantities. The basic flow satisfies 


Ui 

U 2 

P 

P 


u 

oo 

0 

p g 

p g 


0 <_ x 3 <_ » 


(9) 


where all basic flow quantities are assumed constant. The basic flow may be 
regarded as a suitable mean of an actual boundary- layer flow. 


The disturbance equations are as follows: 


i(kgUg - to)u a 
i(kgU 3 - tu)u 3 



dp 

p dx q 

CO o 


i(kgUg - a))p 



i( k gu B ) 



( 10 ) 


j 
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Eliminating u a (a = 1, 2) from the foregoing equations and employing the 
relation 


p = a 2 p 

r oo 

yields an expression for the disturbance pressure in terms of u 3 

_ ip^ (Ug - m/k) (du 3 /dx 3 ) 

P ” ~ [(U g - o)/k) 2 /a 2 ] - 1 


( 11 ) 


where 


u - k 3 U 3 

u g = k 


is the velocity component of the basic flow in the direction of the wave 
normal. When p is eliminated from the resulting equations the disturbance 
equation for the gas is 



k 2 



- m/k ) 2 



OO 


1 


U 3 = ° 


( 12 ) 


Since the real part of w/k equals the phase velocity of the disturbed 
interface, the choice of the proper solutions of equation ( 12 ) depends on 
whether the component of the basic flow in the direction of the wave normal, 
Ug, relative to the phase velocity is greater or less than the speed of 
sound a TO . If it is greater, equation (12) is solved for a disturbance that 
originates at the interface. 

Supersonic solution: 


u 3 = A exp 


ik 


f (Ug - o>/k)' 


- 1 x - 


(13) 


where A is an arbitrary constant. If the phase velocity is less than a^ 
(subsonic), equation ( 12 ) is solved for a disturbance that decays as the 
distance from the wall is increased in the positive direction of X 3 . 

Subsonic solution: 


u 3 = A exp 


(U g - o>/k) 


*3 


(14) 



Ml 
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The disturbance pressure in terms of u 3 for the two cases is: 


Supersonic 


Subsonic 


P = 


p (Ug - 0)/k) 

oo & 

JT ( Ug - u>/k) 2 /a 2 J - 1 


Uc 


ip (U - u/k) 

oo £, 

7 1 - [CUg - ^/k) 2 /a2] 


u 3 


( 15 ) 


(16) 


The kinematic condition at the interface enables one to evaluate the 
arbitrary constant A in the above expressions for u 3 in terms of the ampli- 
tude of the displacement of the interface. Let the equation of the interface 
be 


x 3 (x 1? x 2 , t) = e exp i(kgXg - ait) 
The kinematic condition in the gas at the interface is 


(17) 


u 3 


SXo 9Xn 

TT + U 3 3 ^ 


(18) 


Neglecting second-order disturbance quantities, and evaluating u 3 at 
x 3 = 0 yields 

u 3 = ik (Ug - oi/k)e (19) 


This relation enables one to evaluate the disturbance pressure of the gas at 
the interface x 3 = 0 in terms of e, the amplitude of the displacement. For 
the two cases: 

Supersonic 

ikp^fUg - a)/k) 2 e 

P = ~ — (20) 

7 [ (U g - to/k) 2 /a 2 ] - 1 

Subsonic 

-kp (Ug - w/k)e 

p = ■ (21) 

/l - [(U g - w/k) 2 / a 2 ] 

Subsequently these expressions for p will be used in formulating the 
boundary conditions of the liquid at the interface. For this purpose, it is 
considered sufficiently accurate to neglect the term oi/k compared to the 
term Ug in the above expressions. The expressions that will be used 
subsequently for the gas pressure are: 

8 


L 


Supersonic 


ikpjj|e 
s/M 2 - 1 


( 22 ) 


Subsonic 


~ kp co U g £ 

s/l - M 2 


(23) 


where M h Ug/a^. Note that the Mach number so defined is the projection of 
the actual Mach number of the basic flow in the direction of the wave normal. 
Hence, the flow may be subsonic across the wave crests when the basic flow is 
supersonic. 


Boundary Conditions 

The boundary conditions at the interface are evaluated at the free 
surface. The equation of the interface is 


x 3 ( x i, x 2 > t) = e exp iCk^Xg - tut) 


(24) 


and the kinematic condition in the liquid at the interface is 


u 3 


5x 3 — 9x 3 

3 1 + U P 3Xg 


(25) 


Neglecting second-order disturbance quantities and evaluating u 3 at x 3 = 0 
yields 

u 3 = ik (U^ - iu/k) c (26) 

where is the velocity of the liquid at the undisturbed interface in the 

direction of the wave normal. The first boundary condition to be imposed at 
the interface is that the discontinuity of normal stress equal the surface 
tension. We have for the normal stress 


-P 


g 

1 


3Uo 

+ 2 y 

3x 3 


g 

1 



( 27 ) 


where the second-order terms in the expression for the interface curvature 
have been neglected, and where T denotes the surface tension. 


The above expression evaluated at the disturbed interface is : 



If the effect of viscosity is neglected in the disturbance motion of 
the gas and the relations satisfied by the basic flow are cancelled. 


-P g - Pge + p - 2 u 


8u 3 

3x 3 


Tk 2 c 


(29) 


In the above expression pg is the only term that involves the gas, and all 
quantities are to be evaluated at x 3 = 0. 

Substituting the relations given for p, p g , and e by equations (5), 
(22), and (26), respectively, into equation (29) yields the normal stress 
condition for a supersonic gas flow 


CUi 


co/k) 


du c 

dx. 


u 3 


dU 

dX: 


IV 

+ 

k 


d 3 u 


dx^ 

3 


3 - 3k 2 


dx. 


1 

P 


Tk 2 + pg 


+ lkp g U g N \ 

7m 2 - 1 ^ 


and for a subsonic gas flow from equation (23) 


u 3 

. (U X - co/k) - 

(30) 


du^, 

( U 1 - co/k) 


dU + iv f d3u 3 
dx 3 + k \ dx 3 



The second boundary condition to be 
tangential stress, that is 


= -[ Tk^ + pg 


kp g u g 
si 1 - M 2 


imposed at the interface 


(U x - co/k) 

(31) 

is continuity of 
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I _ f ^ u 3 9u l N \ 

:i3| g = 9 Ti) 


t 23 


g 


f 3Uq 8Up "N 

= + 9^7/ 


l 


( 32 ) 


When the effect of viscosity in the disturbance motion of the gas is neglected. 


Ml i k iu 3 


d 2 ^ 

dx 3 


du n 


£ + 


M ( ik 2 U3 + 


dx c 


du^ 

dx q 


1 = 


= o 


= o 


(33) 


Eliminating u (ct = 1, 2) for the foregoing equation and employing (4) yields 


d 3 


dx 2 

3 


U 3 + k2 _ ke f Vjj , „ 


k Ui - to/k 


(34) 


For a linear velocity profile in the liquid, the last term in equation (34) 
vanishes. In obtaining this boundary condition the following relation 
satisfied by the basic flow has been utilized 


Ul . 

y ir = T 

where t is the shear stress of the basic flow of gas exerted on the 
interface in the direction of the wave normal 


(35) 


k 3 T 63 


Dimensionless Form of the Equations 

The disturbance equations are made dimensionless by a suitable choice of 
reference dimensional quantities. In the present case the depth of the film 
h is chosen as the reference length, and the velocity of the film at the 
undisturbed interface in the direction of the wave normal is chosen as 

the reference velocity; the remaining reference quantity is the liquid 
density p. For a linear velocity profile in the liquid the disturbance 
differential equation (6) is transformed to 
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- 2a 2 4>" + ct 4 (f> = iaR[(F' - c)(4>" - ct 2 4>) ] 


(36) 


where primes denote differentiation with respect to n and where 

F' = U/U} = 1 + n: -1 £ n <. 0 

c = (m/k)/U l 

a = kh 

(j) = (i/kh) (U 3 /UO 

R = U^h/v 

For the supersonic case the boundary conditions at the interface-equations (30) 
are transformed as follows: 


(1 - c)*' “ * - (4> f ' ? - Sa 2 *') = 


1 a z ia 

— + — + 

F 2 W 2 


c.cR/m 2 - lj 


4> 


1 - c 


(37) 


and for the subsonic case 


(l - c)*' - + - ISR - 3 a V) =f -k + ^ 




-r/i - M 2 


1 - c 


(38) 


where 


F = -T 

gh 


a _ 11 


U?ph 


1 = p g U g 

C f R " Pl U2 

F and W are the Froude and Weber numbers, respectively. The boundary 
conditions at the solid surface are 


n = -1; cf> = <J> f = 0 


(39) 
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The differential equation (36) and equations (37) to (39) for the boundary 
conditions constitute an eigenvalue problem for determining the eigenvalue 
c for given values of a, R, F, W, Cf, and M. Another interpretation of the 
parameter c^ can be examined if we note that since x = y(U^/h) 


T = 



(40) 


Numerical Method of Solution of Eigenvalue Problem 

The eigenvalue problem formulated in the previous section is solved in 
the manner employed in reference 7; namely, the parameters a, R, F, W, Cf, 
and M are fixed. Then values of c are found for which equation (36) has 
solutions satisfying the boundary conditions (eqs. (37)- (39)). 

The eigenvalue c and the corresponding eigenfunctions are obtained by 
treating equation (36) as a nonlinear equation. A trial solution is obtained 
by step-by-step numerical integration of the equation starting at n = 0 with 
assumed starting values and with an assumed value of c. The boundary condi- 
tions (eqs. (39), are evaluated at the solid surface n = -1. If the boundary 
conditions are not satisfied, the starting values and c are adjusted by 
means of the Newton-Raphson procedure and another trial is made. To carry out 
the solution of the eigenvalue problem it is convenient to rewrite the fourth- 
order differential equation (36) as a system of second-order differential 
equations. For this purpose, the variable s = <J> !! -a 2 which is related to 
the disturbance vorticity is introduced, and the differential equation (36) 
and the boundary conditions equations (37)- (39) are written in terms of s as 
follows : 

cf) ? r = a 2 cf) + s 

s » t = a 2 s + i a R (1 + n 

at n = 0 

s + 2a 2 cf) = 0 

s' - 2a 2 4> ' - ia^ c z 1 4> - R[(c - l)<t>' + = 0 

at n = - 1 

<|> = <j>' =0 

where for the supersonic case 

_ n ( 1 a 2 A ia 1 

\ F 2 + W 2 / + c f /to 2 - 1 


(42) 


(43) 


(44) 
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and for the subsonic case 


£ = 



a 1 
Cf A - M 2 


( 45 ) 


To carry out the numerical integration of equations (41) four initial condi- 
tions and a value for c must be specified. It is assumed that the other 
parameters have been specified. Since the differential equation and the 
boundary conditions are homogeneous, it is permissible to set <j> = 1 at n = 0. 
This fixes one of the initial conditions. The remaining initial conditions 
may be obtained from the boundary conditions, provided one of them is speci- 
fied. Designate as x the value of <t>'(0). The quantities x and c must 
be adjusted by trial so as to satisfy the boundary conditions at n = -1. The 
initial conditions at n = 0 are given in terms of x and c as follows: 

<t> = 1 

(|> ' = x 
s = -2a 2 

s ' = 2a 2 x + ia^— - j — j - - R[ (c - l)x + 1] 

To satisfy the boundary conditions at n = -1 the values x = <f>'(0) and c are 
adjusted according to the Newton-Raphson procedure. To carry out this proce- 
dure it is necessary to know how <j> and cj>' vary with x and c at n = -1. 

This information can be supplied by formulating the perturbation equations 
(ref. 7) associated with equations (41). The perturbation equations are 
obtained by differentiating the terms in equations (41) with respect to x 
and c. The appropriate initial values for the perturbation differential 
equations are obtained by differentiating the terms in the boundary condi- 
tions (eqs. (46)) with respect to x and c. Integrating the perturbation 
differential equations and the original differential equations (41) enables 
one to evaluate at n = -1 the pertinent quantities: <j>, <f>', c|> x , 4> x ’ , <p c , 

and <p c ’ . The subscripts x and c denote differentiation with respect to 
x and c, respectively. 

The perturbation equations for x are 

<J>x' ' = d 2 <t>x + s x 

s x ' ' = a 2 s x + iaR(F' - c) 
with the initial conditions 

+x = 0 

♦x’ = 1 
S x = o 

s x T = 2a 2 - iaR(c - 1) 
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The perturbation equations for c are 


4>c ? ? = a 2 <j) c + s c 

s c M = a 2 s c + iaR (F ! - c)s c - iaRs 


( 49 ) 


with the initial conditions 


4>c " 0 

*c' = 0 

S r =0 


s c f = -ia 


(C - 1) : 


Rx 


A 


> (SO) 


For values of x and c that yield at n = -1 approximate zeros of t|> 
and <j>*, a better approximation (x + Ax and c + Ac) is obtained by calculating 
the corrections Ax and Ac from the Newton- Raphs on equations 


0 = cp + c£ x Ax + 


(51) 


0 = cj) r + (J)^ r Ax + (f) c ? Ac J 

For given values of the parameters, the procedure outlined should converge to 
an eigenvalue c. 


The results obtained with this numerical procedure will be presented in 
a subsequent section. 


APPROXIMATE ANALYTIC SOLUTIONS 


The region in the neighborhood of R = 0 is of particular interest to 
the problem of wave formation in the melt layer of a body entering the atmos- 
phere since the melt layer is usually highly viscous. Solutions valid in the 
neighborhood of R = 0 may be obtained by an expansion in powers of R as 
suggested by Yih (ref. 8) . However, a somewhat simpler procedure than that 
employed by Yih can be used to obtain solutions for ’’fast 1 ' waves, that is, 
Re(c) >> 1 where Re denotes real part. This procedure consists in neglect- 
ing the variable term n in comparison to the term c - 1 in the differential 
equation (41), solving the resulting equation, formulating the dispersion 
relation in analytic form valid for all values of R, and examining the 
behavior of the dispersion relation in the neighborhood of R = 0. This exam- 
ination will be carried out by expanding the dispersion relation in a Taylor 
series in R. The expansion will be carried out to first order in R. This 
procedure should also yield approximate solutions when |c - l| >> 1. When 
the term n is omitted from equation (41), there are obtained the following 
differential equations with constant coefficients 
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(j> 1 1 = a 2 cf) + s 

s ,! = a 2 s - iaR(c - l)s 


or 


<p r 1 - a 2 (j) = s 

s » * ^ 6 2 S = 0 

where 

B 2 = a 2 - iaR(c - 1) 

the boundary conditions (46) in terms of 3 are as follows: 

Xcj) + (3 2 + a 2 )(J) f - s f = 0 
s + 2a 2 cf) = 0 

where 


la 


c - 1 


f - iaR 


The solutions of equations (52) are 


= Ae an + Be" an + Ce en + De~^ 


s = (B 2 - a 2 )Ce^ n + (B 2 - a 2 )De 


( 52 ) 


(53) 


(54) 


(55) 


(56) 




where A, B, C, and D are arbitrary constants and where B denotes for 
definiteness the root with positive real part of equation (53) . For a non- 
trivial solution of the eigenvalue problem, the substitution of equations (56) 
into equations (43) and (54) requires the vanishing of the following secular 
determinant : 


e -a 

e a 

e-e 


-a 

ae 

-ae 01 

Be~ S 

-Be^ 

X + a(B 2 + a 2 ) 

X - ct(B 2 + a 2 ) 

X + 2ot 2 B 

X - 2a 2 3 

2a 2 

2a 2 

a 2 + B 2 

a 2 + B 2 


After some rearrangement the secular determinant may be written as follows: 
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sinh a 

sinh B 

cosh a 

cosh 8 


a cosh a 

8 cosh 6 

a sinh a 

8 sinh 8 

(57) 

b 2 + a 2 

2a8 

-X/ot 

-X/a 

= 0 

0 

0 

CM 

d 

CM 

a 2 + B 2 



The dispersion relation for R = 0 may be obtained by evaluating the 
secular determinant equation (57) for 3 = a after cancelling the factor 
(3 - ot) 2 from the determinant, and setting 3 = a. The removal of (3 - a) 2 
is required because the secular determinant is satisfied trivially for 3 = a. 
To remove (3 - a) 2 the determinant (57) is rearranged as follows: 


sinh a 

sinh 8 - 

sinh a 

cosh a 

cosh 3 - cosh a 

a cosh a 

8 cosh 8 - 

a cosh a 

a sinh a 

3 sinh 3 - a sinh 

B 2 + a 2 

-(6 - 

a) 2 

-X/a 

0 

0 

0 


2a 2 

6 2 - a 2 


Substituting the identities 

s inh 3 = 

cosh 3 = 


sinh a + 

2 

cosh 


+ 

2 

— sinh 

cosh a + 

2 

sinh 

3_ 

+ 

2 

— sinh 


3 - a 
2 

3 - a 
2 


into the determinant and dividing the second and fourth columns of the 
resulting determinant by the factor (3 - a) yields 


sinh a 


, 3 + a . , 3 - a 

2 cosh — ^ — sinh — ^ — 

3 ^^ 


cosh a 


2 sinh ^ * a sinh 
3 - a 



23 sinh ^ — sinh — ~ — + C osh a 23 cosh sinh - 2 a + sinh a 

a cosh a a sinh a 


3 2 + a 2 


3 - a 
-CB - a) 
0 


-X 

a 

2a 2 


3 - a 

0 

3 + ce 


= o 


(59) 
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letting 0 = a (noting that at 6 = a, 2 sinh[(0 - a)/2]/(0 - a) = 1) yields 
sinh a cosh a cosh a sinh a 

a cosh a a sinh a + cosh a a sinh a a cosh a + sinh a 
2a 2 0 -X/ot 

0 0 2ct 2 


0 

2a 


= 0 


Evaluation of the determinant yields 

(sinh a cosh a - a) (-X/a) = 2a 2 (a 2 + cosh 2 a) 


(60) 


(61) 


With X given by equation (55), equation (61) yields at R = 0 for the sub- 
sonic case 


c ^ _i 1 sinh a cosh a - a 

c_p/l - M 2 a 2 + c °sh 2 a 


(62) 


and for the supersonic case. 


-j 1 1 sinh a cosh a - a 

c - 1 = 2 7~2 (63) 

c^/m 2 1 a + cosh z a 

It should be noted that equations (62) and (63) are valid at R = 0 
regardless of the value of c. This can be seen by carrying out Yih's expan- 
sion procedure for obtaining the dispersion relations. To first order (R = 0) , 
the coefficient of R in the differential equation (41) does not contribute 
to the solution, so that neglecting the term n in comparison to the term 
c - 1 is a superfluous assumption for R = 0. However, this is not the case 
in the neighborhood of R = 0. 

It will be shown subsequently that the dispersion relation in analytic 
form (eq. (59), is extremely useful especially when it is complemented with 
the numerical method previously presented. 

The physical interpretation of the limiting case R = 0 is discussed by 
Yih (ref. 8) in connection with the stability of liquid flow down an inclined 
plane. In the case R = 0 there is no disturbance motion, but merely a sur- 
face corrugation that has no reason to be damped. Since the quantity c is 
expressed in terms of a reference velocity, U^, the actual dimensional phase 
velocity and amplification factor in the case of a surface corrugation would 
be zero since = 0 for a liquid layer with a large but finite viscosity. 

Of considerable significance is the behavior of the dispersion relation 
in this case of static corrugation. Equation (62) indicates that the subsonic 
solutions are unstable at R = 0. This apparently surprising result can be 
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clarified by realizing that the destabilizing action of the gas pressure as 
indicated by equation (23) assumes a dominating role near R = 0 (see eq. (38)). 
Equation (63) indicates that the line R = 0 is a neutral curve for the 
supersonic case. The behavior of c in the neighborhood of R = 0 can be 
determined from equations (53) and (59) . Differentiation of the terms in 
equation (53) regarding R as the independent variable yields at R = 0 


dc _ dc dg . (c - 1) dc 
dR ~ dg dR ~ 1 2 dg 


(64) 


Now equation (59) can be regarded as an implicit equation connecting c and 
g, and the derivative dc/dg can be evaluated at R = 0 by differentiating 
the determinant (59) with respect to 3. This calculation is displayed in 
appendix B. The result for the subsonic case is that the solutions become 
less unstable as R increases from zero, indicating that all terms in equa- 
tion (38) come into play. A further result is that the real part of c 
decreases as R increases from zero, indicating, since Re(c) = 1 at R = 0, 
that the subsonic case does not exhibit the type of instability being examined 
(phase velocity Re(c) greater than unity). Hence, the subsonic case will 
no longer be considered; instead, attention will be devoted to the supersonic 
case only. At R = 0 the result for the supersonic case is as follows: 


dc _ . f (c') 2 
dR 1- ( 4a 

where c ! = c 
Now 


3 _ ^a 2 sinh a (cosh a + a sinh a) 

(a 2 + cosh 2 a) (sinh a cosh a-a) 

(65) 

- 1 and c is given by equation (63). 


c_ 

a 


Cf 


4i 2 -i 


1 

F 2 + W 2 


c (R) = c(0) 



( 66 ) 


and since the right-hand side of equation (65) is purely imaginary, this 
equation gives the damping or amplification of the disturbance in the neigh- 
borhood of R = 0. Setting dc/dR = 0 in equation (65) yields a relation 
among the parameters on the right-hand side valid for neutral disturbances. 
These relations will be presented subsequently, after the approximation 
inherent in equation (65) is evaluated by comparing the results of this equa- 
tion with the numerical results. Before presenting these further results, it 
is of interest to examine solutions of the eigenvalue problem for long waves 
(small a). 

For finite values of R the solution of the eigenvalue problem may be 
obtained in the neighborhood of a = 0 by expansion in powers of a as 
suggested by Yih, reference 8. Solutions of equations (41) are sought near 
a = 0 with c near 1 (eq. (52) indicates that c = 1 at a = 0) . The proce- 
dure set forth for carrying out the numerical solution will be used for this 
analytic solution. The initial conditions to be used are equations (46) . 
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Let 


x = x Q + ax! + . . . 
c = 1 + acj + . . . 

<J> = 4> 0 + + . . . 


5 = 50 + 05 !+. . . 

Substituting these relations in equation (41) yields to lowest order in a 

♦o" - s o 

s 0 " = 0 

the initial conditions equations (46) to lowest order in a are at n = 0 



<j) = 1 


= x r 


s = 


s ' = 


0 

iR 

Ci p2 




( 68 ) 


The satisfaction of the boundary conditions at n = -1: equations (43) 

yield two equations for the determination of the two quantities x Q and 
Ci- This yields the following dispersion relation near a = 0 


c 


ia R 
3 p2 


(69) 


This expression indicates that the line a = 0 is a neutral curve. Further- 
more, it can be verified that dc/dR computed from equation (69) agrees with 
equation (65) to lowest order in a. 


To summarize the results of this section for the supersonic case, we can 
say that the line a = 0 is a neutral curve, that the line R = 0 is a 
neutral curve, that the rate of amplification near the line R = 0 may be 
obtained, approximately, from equation (65) and that the relation between the 
parameters for neutral disturbances near the line R = 0 may be obtained, 
approximately, from equation (65) by setting the left-hand side equal to zero. 
The validity of the approximation inherent in equation (65) will be assessed 
for selected values of the parameters in the following section. 
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Comparison of Approximate Solutions 
With Numerical Solutions 

The numerical procedure described previously was programmed on the IBM 
7094 digital computer at Ames Research Center. Single precision complex arith- 
metic was used. The numerical integration was started with the Runge-Kutta 
formulas and continued with the Adams-Moulton formulas. 


.0002 


The step size was externally controlled and set so that there was no 
significant change in the eigenvalue when the example was rerun with a step 

size equal to half of the original 
value. The step size of 1/32 was 
found to be satisfactory for most runs. 
To assure convergence to an eigenvalue 
the sum of the absolute values of the 
discrepancies of the computed boundary 
values was required to be less than 
5x10” 5 . 



-.0008 


(a) CfA&~- 1 = 0.001 


.018 


.016 k 


.014 
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. 008 ! 


Irn(c) 


.006 


.004 \- 


.002 


-.002 



(b) c f /M 2 - 1 = 0.005 

Figure 2.- Comparison of numerical (dashed line) 
and approximate (solid line) determination of 
a neutral point at R = 0.01 for same value of 
parameters as given in figure 4. 


The approximate analytical 
solutions obtained in the previous sec- 
08 tion were also obtained by means of 
the numerical procedure for selected 
values of the pertinent parameters. 

Because of the large number of 
parameters required to describe the 
physical configuration under study, 
and the ranges of these parameters, 
analytic formulas, such as equations 
(63) and (65), are helpful for describ- 
ing the results of the stability analy- 
sis. Although they may be approximate, 
analytic solutions also simplify the 
searching process usually required to 
find a neutral point. 

The degree of approximation 
inherent in equation (65) may be deter- 
mined from a comparison of the evalua- 
tions obtained from that equation with 
the numerical solutions. These solu- 
tions are compared in figure 2 which 
gives the numerical solutions for 
R = 0.01 and values of Im(c) (Im 
denotes imaginary part) . Also shown 
in figure 2 are the extrapolations 
obtained from equation (66) with 
R = 0.01. It can be seen from this 
figure that the approximate solution 
provides values close to those obtained 
by the numerical method. For example. 
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Figure 3.- Comparison of numerical solutions 
marked (A) with equation (69) shown as a 
straight line at R = 150; Cf»^l 2 - 1 = 0.001, 
U l = 4 cm/sec, h = 0.1 cm, T = 72 dynes/cm, 
g = 980 cm/sec 2 . 


it can be seen from figure 2(b) that 
the numerical solutions yield a value 
of a = 0.187; whereas the approximate 
solutions yield a = 0.185 for a 
neutral disturbance (Im(c) = 0). This 
agreement is considered sufficiently 
accurate for the stability results. 

The parameters used to describe the 
kinetics of the liquid film are the 
depth of the film, h, and the velocity 
of the film at the interface, U^. The 
results shown in figure 2 apply to a 
film with h =0.1 cm and =4 cm/sec. 

The numerical solutions can be 
checked by means of equation (69). 

The results obtained from equation (69) 
are exact for vanishingly small a 
and are compared with numerical 
results at R = 150 in figure 3. 

These two comparisons for 
selected values of the parameters 
lend credence to both the numerical 
and approximate analytic results. The 
results obtained for a more extensive 
range of parameters are presented in 
the next section. 


RESULTS AND DISCUSSION 


The results of a stability analysis are usually displayed in an a - R 
diagram on which the neutral curve (Im(c) = 0) is drawn. As previously stated, 
the line a = 0 is a neutral curve and for the supersonic case, the line 
R = 0 is a neutral curve. However, there are also neutral curves with a 
different from zero emanating from the line R = 0, for example, see figure 2. 
Starting from these points near R = 0, the neutral curves were traced as R 
increased by solving the eigenvalue problem numerically. The results of these 
computations are shown in figure 4. 

An interesting question concerning the neutral curves in figure 4 is why 
does one pass from region of instability to stability as the neutral curve is 
crossed when R is increased at a fixed value of a. This result is contrary 
to the usual situation (Tollmien-Schlichting instabilities). The answer to 
the question lies in clarifying the source of energy for the two types of 
disturbance motion. In the Tollmien-Schlichting instability with c < 1 the 
disturbance energy is supplied by the mean motion of the film. In the 
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present case with c > 1 the 
disturbance energy of the film is 
supplied by the motion of the gas. 

From equation (37) it can be 
seen that as the Reynolds number R 
of the film is increased for fixed 
Cf the influence of the last term 
in the bracket on the right-hand 
side of that equation is reduced. 

The other two terms in the bracket 
always tend to stabilize the distur- 
bance motion for the cases shown in 
figure 4. These two terms become 
dominant when R is sufficiently 
large. 

The last term in the bracket on 
the right-hand side of equation (37) 
is proportional to the pressure coefficient for supersonic flow past a wavy 
wall (ref. 9). Furthermore, the disturbance pressure is not in phase with the 
wall and a wave drag is exerted on the wall (the disturbed gas- liquid inter- 
face in the present case). These phase relations are analogous in a sense to 
the phase relations in the Tol lmien-Schlichting type instability where the 
mechanism for wave generation is the change in phase across the liquid layer 
(ref. 10). However, as opposed to the Tollmien-Schlichting type of instabil- 
ity where a description of the mechanism is rather involved, the mechanism of 
wave generation in the present case can simply be described as supersonic 
wave drag. Because the oblique waves under consideration are three dimen- 
sional, it is not only necessary that the external stream be supersonic, but 
also that the component of the external flow across the crests of the waves be 
supersonic. 

Since these general properties are evident in the a - R diagram for a 
few selected values of the parameters, this representation will not be pursued 
further as a means for displaying the results of the analysis. Rather, 
another representation will be used which is better suited for displaying the 
results for small values of R. 

Because equation (65) provides sufficiently accurate answers for small 
R, the case of interest, it will form the basis for presenting the results for 
neutral disturbances covering a wide range of parameters. Also, it should be 
noted from figure 4 that the neutral curves are nearly horizontal, so that the 
values of a established for neutral points by means of equation (65) near 

R = 0 for a given value of c^/m 2 - 1 should correspond approximately to the 

same value of c^Aa 2 -1 at much larger values of R. 

The values of a corresponding to neutral points may be established if 
one plots equation (65) as a function of a for given values of the param- 
eters and notes where the curve crosses zero. However, it is more convenient 



R 

Figure 4.- Neutral cu rve fo r various values of 
the parameter Cf/M 2 - 1. Region below each 
of the curves is stable, iq = 4 cm/sec, 
h = 0.1 cm, T = 72 dynes/cm, g = 980 cm/sec 2 . 




to set the left-hand side of equation (65) to zero and solve the resulting 

equation. This solution is accomplished by eliminating the term c^/m 2 - 1 
by means of equation (63) and rearranging to obtain the following quadratic 
equation in c f 


a 

2(sinh a cosh a - a) 


3 cosh 2 a + ot 2 


2a 3 cosh 2 


sinh 


a cosn 


- 2 — •) 
a - a/ 



0 


(70) 

As stated previously the kinematics of the liquid film are characterized 
by h and . Given these values and values of a, the surface tension coeffi- 
cient T, and the acceleration of gravity g, all the coefficients in equa- 
tion (70) can be evaluated and the resulting quadratic equation can be solved 
for c f . The choice of the sign before the radical in the quadratic formula 
can be decided by noting that for large a the coefficient of c f2 is 
positive. This fact requires that the positive sign be chosen for positive 
roots (c f positive). After the quadratic equation is solved the parameter 

C£v / M 2 ^ - 1 is evaluated by means of equation (63) . This leads to 


c { /m 2 - 1 = sinh a cosh 

2ac f (a 2 + cosh 2 a) 

The results for neutral disturbances are displayed in dimensional form by plot- 
ting the wavelength of the disturbance A = 27 Th/a against c ^Aa 2 - 1 . These 
results are displayed in figure 5. On each figure the quantity hU-^ is held 
constant. The quantity (l/2)hU^ = Q T is the volume flow rate per unit width 

in the direction of the wave normal. Note that for a given value of cf/M 2 - 1, 
there are either two neutral wavelengths or none. It should be noted that 
the Reynolds number of the film R is not specified for any of the results 
in figure 5. As stated previously these results should hold over a range of 
Reynolds numbers and become increasingly more accurate as the Reynolds number 
approaches zero. For wavelengths much smaller than the depth of the film the 
results in figure 5 may be compared with the results of Willson and Chang 
(ref. 11) who treated the case of an infinitely deep liquid. One of their 
results is that the wavelength is an increasing function of the Mach number. 
This trend is evident from an inspection of figure 5 for short wavelengths. 
However, in contrast to their result, the wavelength does not approach an 
asymptotic value as the Mach number increases. One should be reminded that 
the variables and parameters used in figure 5 correspond to quantities 
measured in the direction of the wave normal, not in the direction of the 
external stream. 
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CONCLUSIONS 


It is concluded from the foregoing linear analysis, that oblique waves 
will be generated in the melt layer of a body entering the atmosphere even 
though the layer is extremely viscous. The mechanism of the wave generation 
is supersonic wave drag. For the type of oblique waves considered, the 
direction of the wave crests is inclined not at the Mach angle of the external 
stream but the direction is such so that there is always supersonic flow 
across their crests. The wavelength of neutral disturbances is governed by 
parameters associated with the external stream and with the melt layer itself. 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif., 94035, Aug. 28, 1968 
124-07-02-23-00-21 
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APPENDIX A 


PRINCIPAL NOMENCLATURE 


a speed of sound 

c phase velocity, dimensionless 

f defined by equation (44) or (45) 

F ! velocity of liquid 

h depth of film 

i imaginary unit 

Im( ) imaginary part 

k wave number 

p pressure 

P pressure, basic flow 

Q f (l/2)U|h quantity of flow per unit width 

Re( ) real part 

s disturbance vorticity, dimensionless 

t time 

T surface tension coefficient 

u velocity component 

U velocity component, basic flow 

x <t>'(0) 

X defined by equation (55) 

Xj,X 2 ,X 3 Cartesian coordinates 


a wave number, dimensionless 

3 defined by equation (53) 

Y ratio of specific heats 



amplitude of disturbance of the interface 
coordinate normal to surface, dimensionless 
wavelength 

dynamic viscosity, liquid 
kinematic viscosity, liquid 
density 

stress components , gas 

disturbance stream function, dimensionless 
angular frequency of disturbance 

Subscripts 

differentiation with respect to c 

gas quantity evaluated at gas- liquid interface 

liquid quantity evaluated at gas- liquid interface 

external flow direction and differentiation with respect to 

tensor index, range (1-2) 

external stream conditions 


total quantity 


Superscript 


Dimensionless Groups 
friction coefficient 
Froude number 
Mach number 
Reynolds number 


Weber number 



K)| h-* 


APPENDIX B 


% 


CALCULATION OF DERIVATIVE 


Differentiation of the determinant (eq. (59)) with respect to 3 and 
evaluation at 3 = a yields 


0 

cosh a 

cosh a 

sinh a 


0 

a sinh a + cosh a 

a sinh a 

a cosh a + sinh a 


2a 

0 

-y i 

0 




a 



0 

0 

2a 2 

2a 


sinh a 

sinh a 

cosh a 

sinh a 


1 a cosh a 

a cosh a + 2 sinh a 

a sinh a 

a cosh a + sinh a 


7 

2a 2 

-2 

-X 

0 




a 



0 

0 

2a 2 

2a 


sinh a 

cosh a 

0 

sinh a 


a cosh a 

a sinh a + cosh a 

0 

a cosh a + sinh a 


2a 2 

0 


0 




dB \a J 



0 

0 

0 

2a 


sinh a 

cosh a 

cosh a 

cosh a j 

1 a cosh a 

a sinh a + cosh a 

a sinh a 

a sinh a + 2 cosh 

a 

2 

2a 2 

0 

-X 

0 



a 


0 

0 

2a 2 

2 
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Evaluation of the four determinants in order yields 


-X 


-4a 4 - 4a 2 cosh 2 a + 2a sinh 2 a 


-4a 3 sinh a cosh a - 2a 2 (1 + a 2 ) cosh 2 a + 4a 4 sinh 2 a 


+2a(sinh a cosh a - a) (- — 


+ (sinh a cosh a - a) ( — 2a 2 (1 + a 


Combining like terms yields 


-4a 2 (2a 2 + 2 cosh 2 a + a sinh a cosh a) 


f X 

+ (— )(2a sinh 2 a + sinh a cosh a - a) 


+ 2a -g- \~ ) (sinh a cosh a - a^ 


SUPERSONIC CASE 


Evaluating equation (55) at R = 0 yields 


X = 


ia 

c - 1 


f 


From equation (44) evaluated at R = 0 


ia 1 
c f /m 2 - 1 


Hence, 


( -X^ _ a 1 _ 2a 2 (a 2 + cosh 2 a) 

\cT/ c - 1 c / M 2 _ \ sinh a cosh a - a 

where the second equality follows after equation (63) is used. 

In order to evaluate (d/dg)X it is convenient to express X 
g rather than R. From equations (44), (53), and (55) we have 


-X = (g 2 - a 2 ) (c - 1) ' 


( 


1 + a 2 > \ a 2 (c - 1) ' 

““"w 2 / cW p-7- 


- (g 2 - a 2 ) (c 


)cosh 2 a = 0 


= 0 (Bl) 


(B2) 


in terms of 

- I )" 1 
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Differentiation of this expression with respect to 3 taking into account 
that c is to be regarded as a function of 6 yields 


1 

a 





a 1 dc 

(c - l) 2 Cf/M 2 - 1 de 

(B3) 


Substituting equations (B2) and (B3) into equation (Bl) and rearranging yields 


(c - 1) dc _ (c 1 ) 2 ^ _ 2a 2 sinh a (cosh a + a sinh a) 

a dg 4a (a 2 + cos h2 ( s i n h a cosh a - a) _ 



This is the desired result, since 


dc . (c - 1) dc 
dR " 1 2 dg 


SUBSONIC CASE 


Substituting equation (45) into equation (55) and evaluating at R = 0 
yields 


-X = ~ia 1 

a C ' 1 c f /l” M 2 


-2a 2 (a 2 + cosh 2 a) 
sinh a cosh a - a 


(B4) 


where the second equality follows after equation (62) is used. 

For the evaluation of the derivative, X is expressed in terms of 3 as 
follows : 

-X = (B 2 - a 2 ) (c - l)- 2 (\ + ^■') + la2 — ■ (c - l)- 1 - (B 2 - a 2 ) (c - l)" 1 

\F Z VT J Cf /l - M 2 


Differentiating this expression with respect to B, regarding c as a func- 
tion of B as before, yields: 


r a i 

(=*X\ 

■ = 2 

O 

r i , 

(± + 

, 1 

ia 1 dc 

dB 1 

l a /_ 

L cc - D 2 ' 

VF 2 W 2 J 

c - 1 

Cf/l - M 2 (c - dP 


(B5) 
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Substituting equations (B4) and (B5) into equation (Bl) and rearranging yields 



The real part of this expression indicates that Re(dc/dR) is negative at 
R = 0. The imaginary part indicates that Im(dc/dR) is also negative at 
R = 0. 
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